Numerical simulation of orbiting black holes 
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We present numerical simulations of binary black hole systems which for the first time last for 
about one orbital period for close but still separate black holes as indicated by the absence of 
a common apparent horizon. An important part of the method is the construction of comoving 
coordinates, in which both the angular and radial motion is minimized through a dynamically 
adjusted shift condition. We use fixed mesh refinement for computational efficiency. 
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One of the fundamental problems of general relativity 
is the two body problem of black holes in a binary orbit. 
Since in general relativity two orbiting bodies emit grav- 
itational waves that carry away energy and momentum 
from the system, the two black holes spiral inward and 
eventually merge. Gravitational waves from black hole 
mergers are expected to be among the primary sources 
for gravitational wave astronomy [l|, |3 ■ 

The last few orbits of a black hole binary fall into the 
strongly dynamic and non-linear regime of general rel- 
ativity, and we therefore turn to numerical simulations 
to solve the full Einstein equations. Numerical relativity 
has seen many advances in recent years, but so far it has 
not been possible to simulate even a single binary black 
hole orbit. The first 3d simulation of a Schwarzschild 
black hole was performed in 1995 In 'i'l, the first 
3d simulation of spinning and moving black holes in a 
"grazing collision" of near-by black holes inside an in- 
nermost stable circular orbit (ISCO) was presented, see 
also . Simulations starting near or even somewhat 
outside an ISCO have been performed e.g. in fRIRI^Hof. 
but after rather short evolution times numerical simula- 
tions of black hole binaries become unstable. In typical 
advanced simulations the evolution time before merger 
is less than 50 Af (where M is the total mass) 0|. An 
open issue is therefore to find methods that allow longer 
lasting evolutions of two black holes before they merge, 
ideally allowing evolution times on the order of one or 
more orbital periods. 

In this paper we present results for a new method to 
choose comoving coordinates that makes it possible to 
evolve two black holes for about one orbital period for 
the first time. The black holes start out close to but well 
outside the ISCO, and the apparent horizons (AHs) do 
not merge before one orbital period has passed. Since 
there are many different choices for the various compo- 
nents of a numerical relativity simulation that crucially 
affect its quality, we will first describe each one of them 
in sufhcient detail to establish our basic framework. We 
will then discuss the major new aspect of our method, 
how we construct comoving coordinates, and discuss our 
numerical results. 

As initial data we choose puncture data 0| for two 
equal mass black holes without spin on a quasi-circular 



orbit based on an approximate helical Killing vector 
[T3 | . Each configuration is determined by the coordi- 
nate distance po of the punctures from the origin. We fo- 
cus on Po = 3.0M, where M is the total ADM (Arnowitt- 
Deser-Misner) mass at the punctures. For po = 3.0M, 
the ADM mass at infinity is 0.985M, the bare mass of 
one puncture is 0.477M, the size of the linear momentum 
of the individual black holes is 0.138M, the angular ve- 
locity is 0.0550/Af, and the orbital period is T = 114M. 
For comparison, post-Newtonian methods and the thin- 
sandwich a ppr oach find the ISCO in the neighborhood of 
T = 65M [iJI, which translates to about po = 1.9M in 
our method. The effective potential method locates the 
ISCO near po = I.IM, T = 35M JJ.,.16J- 

As evolution system we use the modified version of 
the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) sys- 
tem that is described in detail in At the outer 
boundary we impose a radiative boundary condition [TtI 
(we did not implement the monopole term). The black 
holes are handled by introducing a time independent exci- 
sion boundar y ac cording to the "simple excision" method 
described in |l8|. with a generalization from cubical to 
spherical excision regions. We also perform control runs 
without excision using the puncture evolution method 
0, 01 , which typically do not last as long as the excision 
runs, but which allow us to check the excision method. 

As coordinate conditions we use the dynamic gauge 
conditions that proved to be successful for sing le black 
hole runs with and without excision 0, 0, Il9j | and for 
head-on collisions ^Tl\. For the lapse we choose "1-1- log" 
slicing without explicit shift dependence, and for the shift 
we use a particular version of the "Gamma driver" con- 
dition: 
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dtP' = -aP7/'~"BS dtB' = - 7]B\ (2) 

where a is the lapse, is the shift, is its first deriva- 
tive, K is the trace of the extrinsic curvature, is 
the contracted conformal Christoffcl symbol of BSSN, 
and ip is the time independent conformal factor of Brill- 
Lindquist data. After some experimentation we settled 
for our binary runs on m = 4, which helps mimic the sin- 
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gularity avoidance of maximal slicing for puncture runs, 
and for the shift we set n = 2, p = 1, and rj = 2/M. 

One important point to be made about the gauge con- 
ditions C[J-@ is that although they work well for black 
holes without linear momentum, they do not impose 
corotating or comoving coordinates. Moving the black 
hole excision region is showing a lot of promise (29| , but 
here we attempt to minimize the dynamics around black 
holes at fixed coordinate positions by modifying the shift 
condition. Corotatin g fr ames for numerical relativity are 
used for example in [2ll | and with dynamic adjustments 
in 0, |23|. The method that we have developed as a 
first step toward long term comoving coordinates is an 
extension of the methods and ideas of 0, to orbiting 
configurations. 

In order to obtain approximately comoving coordinates 
we introduce the shift vector 

Plom = ^-"{AM-y, ^, 0)* + A2r{-x, ~y, 0)*), (3) 

where a;, y, and z denote Cartesian coordinates, with p = 
(x^ + y'^y^'^ and r = (p^ + z^Y^^. The first term inside 
the brackets is a rotation about the z-axis with angular 
velocity Aiuj, while the second term is an inward radial 
motion with radial velocity A2'rp. The factor ip~'' is used 
to attenuate the shift to zero at each puncture, which is 
needed for simulations without excision. Clearly, for two 
point particles on an inspiraling orbit this shift can cancel 
the dynamics of the point particles completely. For two 
orbiting black holes we can only compensate some aspects 
of the global motion, similar to balancing the bulk motion 
of two stars, with some dynamics remaining in the metric. 

For the runs reported below we have set g = 3, be- 
cause this results in the natural fall off of the shift near 
punctures [l7j . and we use the same value with excision. 
The prefactor Ai can be used to attenuate the angular 
shift for large r ,10] , which simplifies the outer boundary 
and the analysis at large r at the cost of introducing ad- 
ditional differential rotation, but for now we work with 
Ai = 1. Since tends to 1 for r — » oo, the shift 
corresponds to a rigid rotation for large r, in particu- 
lar the coordinate motion becomes superluminal beyond 
a light cylinder. For the radial shift we attenuate with 
A2 = {(? + l)'*/[po(c^ + I Pn)% which is constructed 
such that at the initial radial distance po to the black 
holes the norm of Ai(x,y^^y is unity, at the origin the 
norm is zero, for large p the fall-off is controlled by s, 
and the shape of the attenuation can be adjusted with c. 
We set c = 1 and s = 2. 

To evolve for one orbital time scale it was necessary 
to introduce a dynamic control mechanism with time de- 
pendent velocities and r(i) in the commotion shift 
Q (see also |23)- In order to estimate changes in these 
velocities we define the vector a*(<) = Yl,{^\uncture ~ 
x')Q(t)/ ^ a(t), where the sums run over all points on 
the excision boundary in the orbital plane. The vector 
a*(t) points from the center of the excision region (from 
the puncture) in the direction into which the lapse profile 
has moved off-center. At finite time intervals At, we use 



a*(t) to compute a velocity correction 

Au* = {--ida,npdta\t) - kdr^vea\t)\ At, (4) 

which is designed to damp out motion in a^{t) and to 
drive a^{t) back to zero as in a damped harmonic oscil- 
lator. In coordinates where the punctures are located 
on the y-axis. Aw* defines changes in u}(t) and r{t) by 
Auj = Av^/po and Ar = Av"^. In our case, useful values 
for the coefficients are k^rive = 0.2/M and 'jdamp — 5. 

The evolution of the shift proceeds as follows. We set 
the initial lapse to one and initialize the shift according to 
0, for example with cu — 0.88fi and r = for po — 3M, 
where fl is the angular velocity at infinity defined by the 
initial data. Note that close to the black holes a cor- 
rection to fl is necessary but not unexpected. At each 
time step during the evolution, we evolve the shift with 
0. First, we evolve for a time interval of 5M without 
any commotion correction until lapse and shift have gone 
through their first rapid evolution to adjust themselves 
to the presence of the black holes. After that we compute 
Alli and Ar based on (0J at resolution independent time 
intervals of At = 2M, which defines a shift vector A/3* 
according to This shift vector A/3* is added to /3* 
everywhere on the grid, so the shift changes discontinu- 
ously at intervals of At, but we leave the time derivative 
B' unchanged. 

Assuming a rigidly rotating frame at large distances, 
we generalize the radiative boundary condition taking 
into account that the scalar wave propagation no longer 
happens along the radial direction, and that tensor com- 
ponents have to be rotated to the new frame. For any 
tensor F (indices suppressed) the result is 

d,F^CfiF~v—{F-F^)^^v^—^, (5) 

where C is the Lie derivative, v is the wave speed, and 
Faa is the value of F at infinity. We have experimented 
with cubical and spherical outer boundaries, where the 
latter is expected to have less problems with a global 
rotation. A superluminal shift does not create a prob- 
lem in our runs with the outer boundary at 24M, 48 Af, 
or 96M, since we can lower the Courant factor in the 
outer regions of our fixed mesh refinement grid, which 
we describe below, by switching from Berger-Oliger time 
stepping to uniform time steps. 

All evolutions are carried out with a new version of the 
BAM code ( "bi-functional adaptive mesh" ) |2^ , which is 
built around an oct-tree, cell-centered adaptive mesh ker- 
nel that currently is functional for fixed mesh refinements 
(FMR) without parallelization. Adaptive mesh refine- 
ment (AMR) was made famous in numerical relativity by 
Choptuik's work on critical collapse ,2^, and especially 
in 3d it can offer enormous savings over conventional un- 
igrid codes. However, while the basic technical problem 
of writing AMR codes has been solved many times, see 
e.g. |[25i for an overview and [Hillllllil for some re- 
cent applications in numerical relativity, there have been 
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FIG. 1: Evolution of the AH mass for the black hole binary 
with po = 3.0M. The evolution lasts longer than one orbital 
period of 114M defined by the initial data. The squares mark 
a run with 7 nested levels with coarsest resolution 2M and 
finest resolution h = 0.03125Af , and with the spherical outer 
boundary at about 48M, which crashes around 145M. Also 
plotted are results from seven control runs with the outer 
boundary at 24M and 96M, with a cubical outer boundary, 
and with the AH extracted on a coarser grid to check its con- 
vergence. There is little difference in the results, except that 
the runs with the boundary at 24M last somewhat longer. 



only a handful of examples for the full 3d Einstein equa- 
tions and the evolution of one ^3,[23il23 ^^'^ P black 
holes. The FMR technique with nested boxes of j30j| was 
essential for the feasibility of the first 3d grazing colli- 
sion 0]. 

One aspect of the present paper is that we demonstrate 
that FMR can work successfully even for black holes in 
an orbital configuration. We use nested Cartesian boxes, 
where for black hole binaries with equal mass and no 
spin we only have to store one quadrant of the global do- 
main. BAM's Berger-Oliger FMR algorithm uses third 
order polynomial interpolation in space and second order 
polynomial interpolation in time, following essentially the 
recipe of 0, |30'| . The main missing feature was a reason- 
ably stable unigrid code, which is now available in the 
form of BSSN with dynamic gauge as discussed above. 
An important detail of our setup is the use of the it- 
erative Crank-Nicolson method for time integration. To 
avoid special boundary conditions during Crank-Nicolson 
iterations, BAM uses three buffer points j32i] . 

Let us summarize our numerical results. For the black 
hole binary with po — 3.0M introduced above, evolution 
times of up to 185M are obtained and a typical run easily 
exceeds the orbital period of 114M. Fig.^shows the AH 
mass for one of the black holes as a function of time. It 
is important to note that a common AH enclosing both 
black holes does not form within the achieved evolution 
time, while for sufficiently small values of po (and the 
same AH finder described in j3^ and implemented in 
Cactus [S^l) a common AH is found in [Tol |. 

There is an almost linear drift in the AH mass of about 
10% per lOOM at a resolution of ft. = M/32 near the 
excision region, which becomes smaller with increasing 
resolution as shown in Fig. |21 (We have also evolved 
Schwarzschild on quadrants and full grids for lOOOM and 
more, confirming that our FMR method is convergent in 
the AH mass.) Puncture evolutions without excision give 



FIG. 2: The panel on the left shows convergence of the AH 
mass. The number and size of the refinement levels was not 
changed but the overall resolution was rescaled by a constant 
factor. There is a linear downward drift in the mass which 
becomes smaller with increasing resolution. The panel on the 
right displays the mass at infinity estimated on a sphere of 
radius 20M assuming a Schwarzschild background, showing 
fluctuations of about 20% to 40%. The lower and upper lines 
for a given resolution correspond to a cubical outer boundary 
at 24M and 48M, respectively. 
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FIG. 3: Evolution of the x- and y-components of the shift 
vector along the y-axis. The punctures are located on the 
j/-axis at y = ±3M. 



a quantitatively very similar result, hence the simple ex- 
cision technique does not appear responsible for the drift. 
Since the AH is a slice dependent quantity, the warpage 
of the slice contributes to changes in the AH mass. The 
proper spatial distance between the AHs along the y-axis 
starts at about 9M, rises to IIM, and drops to 7M at 
t — 140M, but since this distance depends on the gauge 
and since it does not converge for the current resolutions, 
this is only a preliminary result. In the future we plan 
to find event horizons to resolve some of the ambiguity. 
Fig.Elalso shows an estimate for the mass at infinity. The 
errors are satisfactory for the present purpose. Since the 
AH mass shown in Fig. ^ is not significantly affected by 
the location of the outer boundary, we conclude that the 
interior of the numerical domain has been computed with 
good accuracy. 

Fig. O shows the evolution of the shift vector. In par- 
ticular, the corotation speed initially increases, then de- 
creases slowly before increasing again towards the end. 
As an indication of the remaining coordinate motion near 
the black holes we show the evolution of the AH in Fig.^] 
A residual drift of similar magnitude is observed also 
for larger separations, which is a likely reason for the 
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FIG. 4: Evolution of the AH of one of the black holes in 
the x-y-plane. The dashed line shows the AH at t — 24M in 
each panel. Initially, the AH moves outward quickly while 
the gauge adjusts itself near the black hole. It then slowly 
shrinks toward the center while being deformed slightly until 
eventually it drifts out of shape before the run fails around 
145M. Note that the proper area changes linearly and only 
on the order of 10% during the entire run, see Fig.0 



code failure that occurs after about 150M rather inde- 
pendently of separation up to po — 12M. 

In conclusion, dynamically adjusted comoving coordi- 
nates enable us to perform the first numerical simulations 
of two black holes near but outside the ISCO for about 
one orbital period. A good indicator for one orbit would 
be the presence of two cycles of gravitational waves. First 
experiments with wave extraction indicate that improve- 
ments of the outer boundary are needed. 
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